rm(list=ls(all=TRUE))

library(foreign)
library(memisc)
library(xtable)
library(Hmisc)
library(maptools)
library(stringr)
library(geosphere)
library(fields)
library(MASS)
library(raster)
library(rgdal)
library(shapefiles)
library(rgeos)
library(multiwayvcov)
library(stargazer)
library(plm)
library(lmtest)
library(data.table)
library(gtools)
library(RColorBrewer)
library(dummies)
library(gtools)


if(Sys.info()["user"]=="IV")
  setwd("/Users/IV/")
if(Sys.info()["user"]=="benjaminbarber")
  setwd("/Users/IV/")


brewer <- brewer.pal(8,"Dark2")
nice_blue <- brewer[3]
fill_blue <- adjustcolor(brewer[3], alpha.f=0.25)
fill_blue2 <- adjustcolor(brewer[3], alpha.f=0.5)
nice_green <- brewer[1]
fill_green <- adjustcolor(brewer[1], alpha.f=0.25)
nice_orange <- brewer[2]
fill_orange <- adjustcolor(brewer[2], alpha.f=0.25)
nice_brown <- brewer[7]
fill_brown <- adjustcolor(brewer[7], alpha.f=0.25)



Combined_Data_Pop <- read.dta("/Users/IV/Dropbox/Nazi/Data/Master Data for Analysis.dta")




#### Creating Raw Distance Plot


Combined_Data_Pop$signal_breaks <- quantcut(Combined_Data_Pop$avg_max_signal_strength_fix, q=20, na.rm=T)
Combined_Data_Pop$signal_quant <- as.numeric(Combined_Data_Pop$signal_breaks)
#Combined_Data_Pop$distance_25_breaks <- quantcut(Combined_Data_Pop$avg_min_dist3, q=25, na.rm=T)
#Combined_Data_Pop$distance_breaks <- as.numeric(Combined_Data_Pop$distance_25_breaks)

pred_signal <- aggregate(Combined_Data_Pop$dec, by=list(Combined_Data_Pop$signal_quant), FUN=mean)
pred_signal <- as.data.frame(pred_signal)
pred_signal[,2] <- pred_signal[,2]*100

#quartzFonts(helvetica=c("Helvetica Light", "Helvetica", "Helvetica Neue Light", "Helvetica Neue Thin"))
pdf("/Users/IV/Dropbox/Nazi/Output/Radio_Singal_Str_Raw_Plot_Decorations.pdf", height=7, width=7)
par(mar = c(5,5,4,1))
plot(pred_signal[,1], pred_signal[,2], ylim=c(15,40), ylab="% of Soliders Recieving a Decoration", xlab="Percentiles of Signal Strength", main="Raw Averages of Decorations Across 
Percentiles of Signal Strength", pch=20, bty="n", axes=F, cex=4, col= nice_blue, cex.lab=1.5, cex.main=1.5)
fixed_labels <- cbind((c(5,10,15,20)), c("25%", "50%", "75%", "100%"))
axis(1, col="white", col.ticks="black", at= as.numeric(fixed_labels[,1]), labels= fixed_labels[,2], cex.axis=1.5)
axis(2, col="white", col.ticks="black", cex.axis=1.5)
#abline(lm(pred_signal[,2]~ pred_signal[,1]), col="darkred", lwd=3)
dev.off()







### Getting Year Dummies for Muster Year ###
myear_dummies <- dummy(Combined_Data_Pop$muster_year, sep="_") #make dummy variables
Combined_Data_Pop <- cbind(Combined_Data_Pop, myear_dummies)




####### Setting Up Variable Lists ######


must_haves <- "+ age_at_muster_relative18 + wounded + nazi2 + urban + NSDAP_Fix + muster_year_1936 + muster_year_1937 + muster_year_1938 + muster_year_1939 + muster_year_1940 + muster_year_1941 + muster_year_1942 + muster_year_1943 + muster_year_1944 + muster_year_1945"
#must_haves <- "+ age_at_muster_relative18 + wounded + nazi2 + urban + NSDAP_Fix + as.factor(muster_year)"
region_fixed_effects <- "+ in_welfare_per1000 + war_per1000 + sozialrentner_per1000 + bigd + logtaxprop + c25juden_share + c25kath_share + c33erlos_share + c33brlos_share + c25anges_share + c25arbei_share"
#population <- "+ population_fix + pop_fix_2nd + pop_fix_3rd + pop_fix_4th + pop_fix_5th"
population <- "+ population_fix"
individual_effects <- "+ catholic + class + married + Relative_Height + Relative_Weight"
sig_str <- "avg_max_signal_strength_fix"
logit_sig_str <- "avg_max_signal_strength_fix_log"
company_fixed_effects <- "+ as.factor(company_id)"
Distance <- "Min_Radio_Dist3"





##### Decorations


decorated_logit <- lm(as.formula(paste("dec", paste(sig_str, must_haves, population, region_fixed_effects, individual_effects, sep=" "), sep="~")), data= Combined_Data_Pop)
summary(decorated_logit)


y_hat <- predict(decorated_logit, type="response")
predicted_values <- cbind(y_hat, decorated_logit$model$avg_max_signal_strength_fix)
predicted_values <- as.data.frame(predicted_values)


radio_range <- seq(min(decorated_logit$model$avg_max_signal_strength_fix, na.rm=TRUE), max(decorated_logit$model$avg_max_signal_strength_fix, na.rm=TRUE), length.out=101)

test_data <- model.frame(decorated_logit)
means <- colMeans(test_data[,3:length(test_data[1,])])


N <- 100000
set.seed(88)
draw <- mvrnorm(N, coef(decorated_logit), vcov(decorated_logit))

predmat <- matrix(NA, nrow=N, ncol=length(radio_range))

for(j in 1:length(radio_range)){
	
	predmat[,j] <- (draw %*% c(1, radio_range[j], means))*100 # multiply them by post-estimated draws

}

lowervec <- apply(predmat,2,quantile,0.025)
uppervec <- apply(predmat,2,quantile,0.975)
medvec <- apply(predmat,2,mean)



adjustcolor("darkblue", alpha.f=0.25)
adjustcolor("darkred", alpha.f=0.25)


pdf("/Users/IV/Dropbox/Nazi/Output/Radio_Singal_Str_Decorated.pdf", height=7, width=7)
par(mar = c(5,5,2,1))
plot(predicted_values$V2, predicted_values$y_hat, type="n", pch=20, col="grey89",  main="Radio Tower Signal Strength & Decoration", xlab="Radio Tower Signal Strength", ylab = "Predicted Probability of Being Decorated", ylim=c(20,55), cex.lab=1.4, cex.main=1.8, cex.axis=1.6, bty="n", axes=FALSE)
axis(1, col="white", col.ticks="black")
axis(2, col="white", col.ticks="black")
#points((radio_range), lowervec, type = "l", lty= "dashed", lwd=1, col="black")
#points((radio_range), uppervec, type = "l", lty= "dashed", lwd=1, col="black")
points((radio_range), medvec, type = "l", lty= "solid", lwd=1, col="darkblue")
polygon(x=c(radio_range, rev(radio_range)), y=c(uppervec, rev(lowervec)), border=NA, col="#00008B40")
#rug(decorated_logit$model$avg_max_signal_strength_fix)
dev.off()


# Less Signal
mean((draw %*% c(1, (mean(decorated_logit$model$avg_max_signal_strength_fix) - sd(decorated_logit$model$avg_max_signal_strength_fix)), means)))

# Signal
mean((draw %*% c(1, mean(decorated_logit$model$avg_max_signal_strength_fix) + sd(decorated_logit$model$avg_max_signal_strength_fix), means)))

# % Change
mean(((draw %*% c(1, mean(decorated_logit$model$avg_max_signal_strength_fix) + sd(decorated_logit$model$avg_max_signal_strength_fix), means)) - (draw %*% c(1, mean(decorated_logit$model$avg_max_signal_strength_fix) -sd(decorated_logit$model$avg_max_signal_strength_fix), means)))/(draw %*% c(1, mean(decorated_logit$model$avg_max_signal_strength_fix) - sd(decorated_logit$model$avg_max_signal_strength_fix), means)))























decorated_dist <- lm(as.formula(paste("dec", paste(Distance, must_haves, population, region_fixed_effects, individual_effects, sep=" "), sep="~")), data= Combined_Data_Pop, family=binomial(link="logit"))
summary(decorated_dist)


y_hat <- predict(decorated_dist, type="response")
predicted_values <- cbind(y_hat, decorated_dist $model$Min_Radio_Dist3)
predicted_values <- as.data.frame(predicted_values)


radio_range <- seq(min(decorated_dist $model$Min_Radio_Dist3, na.rm=TRUE), max(decorated_dist $model$Min_Radio_Dist3, na.rm=TRUE), length.out=101)

test_data <- model.frame(decorated_dist)
means <- colMeans(test_data[,3:length(test_data[1,])])


N <- 100000
set.seed(88)
draw <- mvrnorm(N, coef(decorated_dist), vcov(decorated_dist))

predmat <- matrix(NA, nrow=N, ncol=length(radio_range))

for(j in 1:length(radio_range)){
	
	predmat[,j] <- (draw %*% c(1, radio_range[j], means))*100 # multiply them by post-estimated draws

}

lowervec <- apply(predmat,2,quantile,0.025)
uppervec <- apply(predmat,2,quantile,0.975)
medvec <- apply(predmat,2,mean)




adjustcolor("darkblue", alpha.f=0.25)
adjustcolor("darkred", alpha.f=0.25)


pdf("/Users/IV/Dropbox/Nazi/Output/Radio_Min_Dist_Decorated.pdf", height=7, width=7)
plot((radio_range), medvec, type="n", pch=20, col="grey89",  main="Radio Tower Distance & Decoration", xlab="Distance from Closest Radio Tower (in km)", ylab = "Predicted Probability of Being Decorated", ylim=c(20,40), cex.lab=1.4, cex.main=1.8, cex.axis=1.6, bty="n", axes=FALSE)
axis(1, col="white", col.ticks="black")
axis(2, col="white", col.ticks="black")
points((radio_range), medvec, type = "l", lty= "solid", lwd=1, col="darkblue")
polygon(x=c(radio_range, rev(radio_range)), y=c(uppervec, rev(lowervec)), border=NA, col="#00008B40")
#rug(decorated_logit$model$Min_Radio_Dist3)
dev.off()



# Close
mean((draw %*% c(1, 0, means)))

# Far
mean((draw %*% c(1, 100, means)))

# % Change
mean(((draw %*% c(1, 100, means)) - (draw %*% c(1, 0, means)))/(draw %*% c(1, 0, means)))


